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ABSTRACT 

Aims. Accurate measurement of the cosmic microwave background (CMB) anisotropy requires precise knowledge of the 

instrument beam. We explore how well the Planck beams will be determined from observations of planets, developing 

techniques that are also appropriate for other experiments. 

Methods. We simulate planet observations with a Planck-\ike scanning strategy, telescope beams, noise, and detector 

properties. Then we employ both parametric and non-parametric techniques, reconstructing beams directly from the 

time-ordered data. With a faithful parameterization of the beam shape, we can constrain certain detector properties, 

such as the time constants of the detectors, to high precision. Alternatively, we decompose the beam using an orthogonal 

basis. For both techniques, we characterize the errors in the beam reconstruction with Monte Carlo realizations. For a 

simplified scanning strategy, we study the impact on estimation of the CMB power spectrum. Finally, we explore the 

consequences for measuring cosmological parameters, focusing on the spectral index of primordial scalar perturbations, 

ris. 

Results. The quality of the power spectrum measurement will be significantly influenced by the optical modeling of the 

telescope. In our most conservative case, using no information about the optics except the measurement of planets, we 

find that a single transit of Jupiter across the focal plane will measure the beam window functions to better than 0.3% 

for the channels at 100-217 GHz that are the most sensitive to the CMB. Constraining the beam with optical modeling 

can lead to much higher quality reconstruction. 

Conclusions. Depending on the optical modeling, the beam errors may be a significant contribution to the measurement 

systematics for Ug. 

Key words. Cosmology: cosmic microwave background - Cosmology: cosmological parameters - Cosmology: observations 
- Data Analysis: methods 
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\^ ' 1. Introduction nal of the covariance matrix. Because of the conservative 

t:;j- treatment of beam errors in the three-year release likeli- 

CT) ' Robust measurements of the cosmic microwave background hood method, the cosmological estimates fortunately do 

j^ ; (CMB) anisotropy, the source of much of our understanding not change much — mostly manifesting as a 0. 7cr shift in the 

^ . of the universe's contents, geometry, and primordial fluctu- present-day amplitude of pertu rbati ons, a^ ([Spergel et all 

(^ ations, require detailed control over the systematics of the |2007l : iDunklev et all 120091 : iKomatsu et al.ll2009[ l 

,-H instrument. The spatial response to a signal on the sky. In this wor k, we examine the recently launched 

r^ known as the point-spread-function (PSF) or simply the Plan ck missior£l (iTauber et al.ll2010l : iPlanck Collaboration! 

.1^ telescope beam, is an important systematic effect because |2006^ . the next generation satellite to measure the CMB 

K> it smooths the anisotropy on the sky, damping high spa- anisotropy. To wring the full cosmological information from 

^ ■ tial frequencies in the angular power spectrum and washing the observations, the proper calibration of the beam over a 

C3 ■ out the encoded cosmological information. To recover the wide range of angular scales will prove crucial. Because the 

power spectrum, we face the challenging task of accurate sensitivity of the detectors alone would allow a cosmic vari- 

beam reconstruction. ance limited measurement of the temperature power spec- 

The release of the five-year results from the Wilkinson ^rum to high multipoles, determination of the physics at 

Microwave Anisotropy Probe fWMAP) (iHinshaw et ap high spatial frequency will depend on the removal of sys- 

2009) highlights the issue's importance. iHill et alJ 0oM tematics, in particular the quality of the beam reconstruc- 



substantially refines the model of the instrument beam over tion, especially since errors in the beam imprint errors on 

the previous versio n, which is then f olded into the power ~; ~~ ; " "" . /r^, , x . . r i 

, ^- ^ Jitm 7 — nionnnh r^^ i^ • • Planck [nttp: / /www.esa.int/Planck) is a proiect ot the 

spectrum estimate (INolta et al.l 120091 ). The result is an m- t^ o a t:^o a -^i • ^ ^ -j j u 

^ -.in ^ — \ — .1 .1 buropean bpace Agency — bbA — with instruments provided by 

crease m the five-year power spectrum over the three-year ^^^ ^^i^r^ii^^ Consortia funded by ES A member states (in par- 

of 2 percent at the first acoustic peak and slightly more ^icular the lead countries: France and Italy) with contributions 

at smaller scales. These changes are easily visible by eye from NASA (USA), and telescope reflectors provided in a col- 

when plotting the three year and five year spectra together, laboration between ESA and a scientific Consortium led and 

and outside the nominal error bars taken from the diago- funded by Denmark. 
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the power spectrum which are strongly correlated between 
multipoles. 

The beam error thus will affect a diverse range of 
science goals for the Planck CMB maps and power 
spectra. These include constraints on the early uni- 
verse, in the measurement of the primordial spec- 
trum's slope and running, the CMB da mping tail, or 
any exotic physics at recombination (e.g [Colombo et alJ 
l2QQ9l ). In the later universe the high-/ spectrum af- 
fects constrain ts on the ma tter distribution from CMB 
lensing fe.g. Amblard et al. 2004: Kesden et al. 20031: 
i Hirata fc Seliakl l2003bl : iLewisI [20051 : iLewis fc Challinorl 
12006") , lensi ng-derived limits on neutrino masses from 
CMB alone (jKaplinghat et al.l 120031 : llchikawa et aD I2005D 
and i n combination with o ther larg e scale st ructure 
data (iKristiansen et al.l 120061 : ide Bernardis et al.l [20081 
information on cluster physi cs from the Sunyaev- 
Zeldovich (SZ) power spectrum ("Komatsu fc SeljaH [20021 : 
[Diego fc Maium dar 2004; Holder et al. 2007), and models 
of c orrelations in point source populations (e.g. [Righi et al.[ 
[2008) . Finally, uncertainty in the beam shape adds error 
to cluster SZ and point-source flux measurements. In addi- 
tion, understanding Planck^ s beam error is important for 
other experiments when forecasting cosmological perfor- 
mance based on Planck prior parameter constraints. 

Historically, CMB experiments have used a combination 
of optics calculations and plan et measurements t o work out 
the shape of the beam (e.g. [Page et aD [2003[ : [Crill et aD 
[2003[ : [Masi et al.[[2006l among many others). Planets prove 
so useful because as bright, compact sources, they resemble 
^-function signal impulses. In Sec. [2l we discuss Planck^ s 
planet observations during the course of routine opera- 
tions, our pipeline for simulating planet observations, and 
two methods for measuring the structure of the instrument 
beam: one in which we use significant prior information 
about the beam's shape, and another where we use very lit- 
tle. We then interpret these beam reconstructions in terms 
of their effect on the CMB power spectrum. In Sec. [3l we 
present the results of our computations, including detailed 
forecasts for the characterization of the beam reconstruc- 
tion and uncertainties. In Sec. [H we examine the impact 
on the scalar perturbation spectral index, n^. Finally, we 
conclude in Sec. [5l 



2. Methods 

2.1. Scan strategy and planet properties 



The Low-Frequency Instrument (LFI, [Bersanelli et al.[ 
2010[) and t he H igh- Frequency Instrument (HFI, 
Lamarre et al.[ [20101 ) of the Planck spacecraft will 
observe the CMB from an orbit around the Earth- Sun 
Lagrange point, L2. The spacecraft spins at ~ 1 rpm with 
the spin axis pointed roughly in the anti-Sun direction, 
sweeping beams (oriented 85° from the spin axis) across 
the sky. The spin axis is stepped along the ecliptic, 
about hourly, to keep the spin axis pointed away froni 
the Sun. Th e detailed strategy ([Dupac fc Tauberl [2005[ : 
[Tauber et al.[ [2Q10) modulates the spin-axis direction in a 
cycloid pattern, yielding virtually complete sky coverage 
in 7.5 months. 
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Fig. 1. Simulated Jupiter observation with a Planck 100 
GHz horn. Each point represent a single sample of the time- 
ordered data, and the coordinates are in the planet frame, 
with the X-axis in the cross-scan direction. In the left panel, 
the telescope pointing at each sample is clearer, and the 
half- maximum curve is marked. 



Since the anti-Sun direction and planet ephemerides 
may be predicted well in advancqj, it is straightforward 
to estimate the times for planet observations by Planck. 
Because the exact observation time of a planet by an indi- 
vidual detector depends on the detector's position within 
the focal plane, the details of Planck^ s orbit around L2, and 
constraints on the spin axis modulation, our estimates are 
correct to ~ 1 week. Planet brightness is determined by the 
orbital configuration during Planck observation. Since the 
outer planets orbit the Sun more slowly than Planck, they 
will be observed roughly once per sky survey. The orbital 
geometry required for observation (the planet is seen from 
L2 at ±85° from the anti-Sun ray) means the distance to 
the planet, and hence the brightness, is similar at every 
observation. The planets will be at high galactic latitude 
during their observation by Planck in 2009 and 2010, so we 
have not included galactic emission here. 

The data from a single observation of a planet will con- 
sist of a consecutive set of measurements for the duration 
of the focal plane's transit over the planet. Since planets lie 
near the ecliptic, the detector pointings in a single observa- 
tion fall in stripes which are approximately perpendicular 
to the ecliptic. 

Here we consider two realistic effects on the pointing 
of th e spacecraft which reflect dynamics simulations from 
ESA (Tauber et al.''2010). First, every re-pointing will differ 
from the desired pointing by a small random error. Second, 
the spacecraft spin axis nutates with a ~ 6 minute period 
and an amplitude that may change after each re-pointing 
of the spacecraft. The nutation of the spacecraft spin axis 
spreads samples along in cross-scan direction. We based 
these effects on pointing simulations, but final values will 
be determined in flight. Figure [T] shows a simulation of 
a Jupiter observation with P/ancA:-like pointing, using the 
pipeline developed for this work. 

Estimates of the brightness temperatures for each 
planet in each band are shown in Table [H Using angular 
diameters from the ephemerides, brightness, and assum- 
ing nominal beam sizes, we can calculate the beam dilu- 
tion factors to find the signal seen by each channel. Jupiter 
and Saturn will be extremely bright in all the bands. Mars, 



^ Available online via the JPL HORIZONS system: 
http : //ssd . jpl . nasa . gov/?horizons 
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Uranus and Neptune will be seen with high signal to noise 
and will be useful for the main lobe of the beams (as may 
thousands of galactic and extra-galactic compact sources). 

2.2. Beam fitting methods 

We have experimented with several methods to fit the 
beams. Our original efforts concentrated on fitting beams to 
simulated, noisy maps of planet observations, a procedure 
which we ultimately found unsatisfactory. We made maps 
two ways: binning on a rectangular grid after an offset sub- 
tractio n (to remove long-time drifts) and onto a HEALPix 
griqj ^Gorski e t al."2QQ 5D using the destripin^ mapmaker 
Springtide (e.g. Ashdow n et al.ll2009[ ), and attempted to fit 
beam parameters by minimizing x^ over the pixels. Initial 
fits gave rough beam parameters, but they were incorrect 
in detail due to the pixelization effects. For the smaller 
beams, map pixels near bright planets contain large signal 
gradients (at the resolution of the Planck CMB maps). In 
principle increasing the map resolution solves this problem, 
but the number of pixels required approaches the number 
of time samples. Therefore we abandoned map-domain fit- 
ting for time-domain fitting, the focus of this current effort 
(and also employed successfully by Burigana et al. 2000). 
This approach also proves convenient for studying some 
effects (e.g. noise correlations, bolometer time constants) 
which are more easily represented in the time or frequency 
domain. 

In passin g, we mention a c lever alternative method 
proposed by IChiang et al.l (|2002[ ) , where the asymmetric- 
beam-induced statistical anisotropy of the observed noisy 
CMB field is compared in Fourier space to statistically 
isotropic noise realizations to deduce the beam asymmetry, 
but not the complete beam window function. 



the planet on the sky, appropriate for the few-hour time 
scales h ere. 

Taub er et al.l (|2010f ) summarizes the optical properties 
of the Planck mission's telescope. For the beam, we use the 
detailed calculations produced by the collaboration based 
on models of the telescope optics (ISandri et al.ll2QQ2l I2OIOI : 
Yurchenko et al. 2004:M affei et aDl2010[ ). Beam values are 
provided on a tabulated grid, which we evaluate at non-grid 
points by interpolation. We fit a Gaussian to capture the 
beam's largest scales, then use 2-d cubic spline interpolation 
to reproduce the residuals to this fit. The interpolated beam 
is the sum of the Gaussian and the spline interpolation and 
reproduces the gridded beam exactly on pixel centers. We 
model the planet as a point source, so that after convolution 
with the beam, the planet signal resembles the beam shape, 
with peak temperature given by Table [H 

To include the impact of the CMB on the planet fits, 
we simulate small scale CMB modes. These are computed 
by FFT in a flat sky approximation on a plane surrounding 
the planet scan, expanded to avoid edge effects in the beam 
data. Because of the high planet signal, we find that the 
CMB does not have a material effect on the beam recovery. 



2.4. Detector properties 

Our simulated detectors are primarily characterized by 
their noise attributes, which we set to mimic the actual 
Planck detectors (see Table [2]). Optionally, for the HFI de- 
tectors, we include a time-constant and/or nonlinear re- 
sponse in our simulations. 

To capture low-frequency drifts in the electronic am- 
plifiers and bolometer temperatures, we use a noise power 
spectrum of the form 



P„(/) = Pwhite[l + (///k„ee)-"] 



(1) 



2.3. Rapid IVIonte Carlo simulation 

We designed and implemented a software pipeline to rapidly 
simulate planet crossing in the time domain, and then re- 
construct the beam. The Planck Collaboration has im- 
plemented an extensive s oftware infrastructur e to simu- 
late time-ordered data (e.g. lReinecke et ani2006[ ), but these 
tools are, by design and optimization, intended to simulate 
large surveys. Here we want to examine the beam fitting 
procedure by the Monte Carlo method, focusing our inter- 
est, by contrast, on the small fraction of the data near the 
planets, which permits optimizations in the design of our 
pipeline suited to that task. On a laptop, our pipeline can 
simulate a planet crossing and subsequent beam reconstruc- 
tion in a few seconds, fast enough that, on a cluster, we can 
rapidly generate thousands of simulations. The modeling 
includes simulated pointing, realistic beams, planets, 1// 
and white noise, the CMB, and several time-domain filters. 
We test our beam fitting methods using these simulated 
observations, characterizing the beam errors by the Monte 
Carlo method. 

The pointing is generated on rings and includes a ran- 
domized re-pointing error between rings. We model nuta- IPel abrouillel 1 19981 
tion of the satellite spin axis as a cross-scan oscillation at 
a fixed frequency. We translate the pointing into the frame 
where the planet is fixed, accounting for linear motions of 



where the low frequency index a ~ 1.7 for LFI and a ~ 2.0 
for HFI. We consider this noise as a sum of correlated and 
white parts, generated in separate steps. The correlated 
low frequency part is generated via an FFT, and is con- 
tinuously but slowly sampled (typically ~ 1 Hz) for the 
duration of the planet crossing (16-96 hours, depending on 
the beam size), then interpolated to the detector sampling 
frequency (up to 200 Hz) only when the detector is close to 
the planet. The white noise is sampled at the detector rate, 
but generated only near the planet. This multi-scale ap- 
proach is much faster than generating the noise at the full 
data rate for the duration of the crossing. The interpolation 
and slow sampling of the correlated noise realization causes 
a smoothing of the low frequency portion of the noise, but 
the slow sampling rate is chosen based on the knee fre- 
quency so that the white noise masks this smoothing, and 
yields Gaussian noise with a very close approximation to 
the desired power spectrum. 

In practice, Planck^ s data streams will be filtered to 
decrease the impact of the low frequency noise, particu- 
larly in the course of mapmaking. One promising way to 



achieve this is th r ough destriping JB urigana et al.l 



Kurki-Suonio et al. 



Revenu et al.l 1 1998 : Sbarra et al.l 



1997 : 



20031 : 



12004 ). which involves fitting offsets to 



http : //healpix . j pi . nasa . gov 



the noise, using crossing points in the sca n as points ofref- 
erence to separate signal and noise. Te renzi et al.l (|2004l ) 
found for LFI detectors that undestriped 1// imparts 20- 
30 percent systematics to the fiux recovery of 1 Jy point 
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Band 
(GHz) 



Jupiter^ Saturn^ 



Peak temperature 



Mars 



(mK) 
Uranus^ 



Neptune^ 



30 


44.8 


7.72 


2.41 


0.290 


0.194 


44 


90.6 


15.3 


4.75 


0.564 


0.319 


70 


303 


49.3 


15.0 


1.78 


0.830 


100 


754 


121 


38.1 


4.29 


1.73 


143 


1.86 X 10^ 


299 


91.8 


9.12 


3.66 


217 


7.58 X 10^ 


1.22 X 10^ 


350 


31.9 


12.7 


353 


3.59 X 10^ 


5.60 X 10^ 


1.67 X 10^ 


131 


51.9 


545 


3.72 X 10^ 


6.50 X 10^ 


2.31 X 10^ 


1.59 X 10^ 


630 


857 


3.31 X 10^ 


5.14 X 10^ 


1.89 X 10^ 


1.11 X 10^ 


4.43 X 10 



Table 1. Model pea k te mperature (CMB t h ermody n amic unit s) in the Planck beam for planets. References: (1) 
iNaselskv et J] (|2QQ7D : (2) iGoldin et all (|1997D : IWrightl (|1976[ ): (3) iGriffin fc Orto] (|1993D . 



sources (roughly Neptune's flux at 44 GHz). The noise af- 
ter destriping may be characterized with an effective powe r 
spectrum. Using the analysis of Kurki-Suonio et aU (|2009[ ), 
we express the power spectrum of this noise as 



Pd. 



estripe 



(/) = ^n(/)x 



10" 



sin^(7r/toff) 
(7r/toff)2 



(2) 



where the offset duration parameter, toff, produces the low- 
est noise residuals near l/(2/knee) (the precise minimum of 
the residuals depends somewhat on the details of the signal 
and scan). We consider this noise as the residual 1// in our 

planet observations. 

The HFI aboar d Planck (Lamarre et al.ll2003[ ) consists 
of 52 bolometers ([Holmes et al.l |2008[ ) fed by feed horn 
structures, read out at nearly 200 Hz. The bolometer's ther- 
mal response to an the incoming optical signal is described 
by a transfer function expressed in the Fourier domain as 
a single pole low-pass filter: 



T(u;) 



1 



1 



(3) 



where cu is the angular frequency of the signal and r is the 
thermal time constant of the bolometer. In analysis of CMB 
data, the detector time constant can be treated as part of 
the effective beam or simply deconv olved from the tim e 
ordered data as a pre-processing step ('H ananv et al.lfl998[ ). 
Here we treat the bolometer time constant as an additional 
parameter in our model of the instrumental response to 
planet observations. 

In practice, the details of the bolometer's thermal circuit 
can lead to a transfer function that is not described by 
a single-pole low-pass filter. In principle we can include a 
more general transfer function as additional parameters in 
our fit. 

The ambient optical background of the HFI bolometers 
is dominated by the CMB and thermal emission from the 
telescope and optical filters. The planets Mars, Jupiter, and 
Saturn are expected to be extremely bright compared to 
the ambient optical background and will drive the bolome- 
ters nonlinear. In th e case of HFI's readout electronics 
(jGaertner et al.lll997[ ), a drop in bolometer resistance al- 
ways leads to a state of "overcompensation" for the tran- 
sient compensation. The average signal over a readout cycle 
will not drop as much as in a DC biasing scheme, mitigating 



the effect of bolometer nonlinearity. To simulate the nonlin- 
ear response of HFI, we use a gain curve fro m a model of the 
detec tor and readout electronics (DESIRE, ICatalano et al.l 
I2009D . 

In practice, the response can be modeled and corrected, 
however when we include nonlinearity here, we will simply 
cut the samples most affected by nonlinearity, which should 
give a performance baseline that we should exceed. 



2.5. Beam model I: parametric linear distortions 

In the final analysis of Planck data, the observations of 
the planets can be used to constrain the principal com- 
ponents of a parameterized beam model, based on optical 
computation of the system of mirrors and horns. At suffi- 
cient accuracy, this model defines a family of beams which 
can faithfully represent the actual beam. Though the full 
analysis is beyond the scope of this work, we can proceed 
fruitfully using a linear approximation to the beam distor- 
tion, which can be described by seven parametersQ 

When tracing rays, a distortion in the optics can be 
represented by a transformation of the ray destinations. A 
ray from the source which originally arrived at the image 
plane at x will arrive at x' = T(x) after distortion. For 
nearby rays and small distortions, we can expand in a series 
to linear order: 



X 



T(x) 



Tix 



(4) 



where Tq is a vector in the plane and Ti is a 2 x 2 matrix. 
This is equivalently expressed as 



x'«Ti(x-xo) 



(5) 



where we introduce xq as a beam offset in preference to Tq. 
Therefore, if 5true(x) defines the realistic beam at po- 
sition X, then the distorted beam can be written as 



^model(x) = A ^true (x ) 



(6) 



where A is the relative amplitude. It is more convenient to 
work with the inverse of Ti and decompose its four ele- 
ments into a rescaling, a rotation, and two components of 



^ Very simple parametric models of the beam, like an ellipti- 
cal Gaussian, cannot faithfully represent the beam simulations, 
and impart noticeable biases to the beam window function. See 
Sec. [2:51 
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Band Sample rate White noise Knee for 1// Low-/ Bolometer FWHM 

(GHz) (Hz) (/iK s""*^^^) (mHz) index a r (ms) (arcmin.) 

Low-Frequency Instrument 



30 


32.5 


170 50 


1.7 




32 


44 


46.5 


200 50 


1.7 




20 


70 


78.8 


270 50 

High- Frequency Instrument 


1.7 




13 


100 


185 


50 30 


2 


10.3 


9.2 


143 


185 


62 30 


2 


4.5 


6.5 


217 


185 


91 30 


2 


3.2 


4.5 


353 


185 


277 30 


2 


4.2 


4.2 


545 


185 


1998 30 


2 


1.5 


4.2 


857 


185 


91000 30 


2 


1.9 


4.0 



Table 2. Planck detector properties used for modeling. The HFI sampling rate is approximate and subject to on-orbit 
tuning. Noise figures give values in C MB thermo dynamic units. F WHM based on a Gau ss ian fit to the realistic beam 
models mentioned in text. References: iPlanck Collaborati on (2 0061) ; lAshdown et al.l (|2009f ): iHolmes et al.l (|2008.) . 



shear: 



1 + s 

1 + 5 



x(l 



■ll 



R(V^) 

2 \-l 



1 + 7+ 7x 
7x 1-7+ 



(7) 



■7 



where s defines the rescahng, R(V^) is a rotation through 
angle ip^ and 7+ and 7x are respectively the perpendicular 
and diagonal components of shear. 

Note that this simple parametric family contains the 
realistic beam, when the offset xq is zero and the transfor- 
mation matrix Ti is the identity. Because the rotation and 
shear transformations have unit determinant, the s param- 
eter completely characterizes the solid angle of the beam, 
with ft oc {1 -\- s)'^ . Because of the comparative rigidity of 
this model, we find that we can successfully constrain de- 
tector transfer function parameters or multiple planet am- 
plitudes as additional parameters in a single fit. 

This parameterized model, with as little as seven pa- 
rameters per beam, is appeal ing, but probably be too sim- 
plistic. For example, in the iHill et ahl (|2009[ ) analysis of 
WMAP^ r^ 430 parameters are used to fit simultaneously 
the ten beams of each telescope. Our parametric model's 
simplicity implies a rigidity, which manifests itself in prob- 
ably too optimistic beam errors. Later we compare with a 
more flexible model. 



2.6. Model parameter fitting and convergence 

To fit the parameters to the data samples d^, we minimize 



X 



7 , [di — 5model(Xi)] C^j [dj — 5model(Xj 



(summed over time samples i, j) with a downhill simplex 
method. There are a few 10^ time samples within a few 
beam full- width half-maximums (FWHMs) of a planet. 
This makes fitting in the time domain tractable for a sparse 
covariance matrix, but possibly not for a dense matrix. We 
only considered diagonal matrices, appropriate for white 
noise. This is a simplification when we include 1// noise 



and CMB. In practice this does not to bias our result (av- 
eraged over an ensemble of noise and CMB realizations), 
although it makes the errors on the fit larger than with an 
optimal estimate. 

In our parametric beam model, the vector of parameters 
which exactly recovers the true beam is 



{A, Xo, S, V^, 7+, 7x } = {Apianet, (O', O'), 0, 0°, 0, 0}. 



(8) 



Because of noise and sparse pointing, the set of parameters 
which minimizes y^ will differ from these, and the distri- 
bution of these parameters characterizes our uncertainty 
in the beam recovery. The sources of noise (detector and 
CMB) are Gaussian, so the y^ minimization will seek an 
unbiased estimate for the beam function at each sample. 
The beam parameters are related to the beam via a non- 
linear function, so they can be biased in our estimate, even 
when the beam is unbiased. In practice, however, we find 
these biases to be negligible. For example, fitting to Jupiter 
after '^10^ MC steps for our highest planet-signal-to- noise 
channel (857 GHz), the mean cross-scan position is about 
1.5 times the predicted standard deviation of the mean for 
this sample size. This bias is 5 percent of the typical disper- 
sion for a single realization, and only 10~^ times the beam 
FWHM. Other channels show similar biases, and the biases 
on other parameters are similar or better. 

At each step in the Monte Carlo simulation (each with 
an independent realization of pointing, noise, and CMB), 
we start the minimization at a random point in parameter 
space, located with uniform probability within a rectangu- 
lar solid centered on the parameters of the true beam, with 
dimensions {0.5 x Apianet, (2', 2'), 0.2, 10°, 0.2, 0.2}. The di- 
mensions of this box, cut in half, are used for the step size 
to initialize the simplex minimizatioE[j. We execute the al- 
gorithm for a fixed number of iterations, which sets the 
accuracy of the minimization at each step in the Monte 
Carlo. Parameter distribution convergence depends on this 
accuracy and on the total number of Monte Carlo steps. 



^ We use an implementation (nms implex) from the GNU 
Scientific Library, version 1.10 
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To verify the convergence of our beam parameter distri- 
butions, we ran a suite of simulations, varying the iteration 
count in the minimization. We consider one horn from ev- 
ery band of Planck observing Jupiter, with destriped noise 
and CMB, including a time constant for the HFI detectors. 
Depending on the channel, we find that after 1,000-3,000 
iterations, the simplex minimization will have settled suffi- 
ciently to not have a material effect on the variance of the 
parameter distributions. (More parameters require more it- 
erations to converge.) The 353 GHz beam is the slowest to 
converge in the simplex fitting, and the only channel which 
occasionally settles in spurious local minima in x^ (more 
fully discussed in Sec. 13. ip . The variance of beam parame- 
ters converges after roughly 1000 Monte Carlo steps. 

2. 7. Beam model II: non-parametric basis functions 

We have argued that our parametric model may be too rigid 
to give realistic beam errors. As an alternative, we want to 
constrain the beam directly by the planet measurements, 
without recourse to an optical model. Parametric methods 
are more powerful statistically, but non-parametric recon- 
struction should provide a robust consistency check. In ad- 
dition, if the beams contain tails, corrugations, or other fea- 
tures on-orbit which were not present in ground measure- 
ments or simulations, the non-parametric method provides 
a way to represent them. 

To represent an arbitrary beam, any complete basis will 
suffice, but some represent the beam more compactly. The 
simulated beams of most channels are approximate ellipti- 
cal Gaussians, so the eigenfunctions for an asymmetric 2-D 
quantum harmonic oscillator form a convenient basis. We 
construct the ground state of the Hermite-Gauss functions, 
an elliptical Gaussian, so that it minimizes the square de- 
viation from the simulated beam. Below we illustrate the 
computation of basis coefficients with integration on the 
sparsely sampled sky using the orthogonality relation, as 
opposed to fitting for parameters. 

Similar basis functions (but based on an axisymmet- 
ric Gaussian) are frequently used in astronom y, for ex- 
ampl e to describe galaxies shap es ("shapelets," iRefregierl 
120031 : iMassev fc Refregierl |2QQ5|) and telescope PSFs in 
gravitational lensing studies (Bernstein & Jarvis' 2002; 
iHir ata & S eliak 2003a) . In particular, Refregier (2003) con- 
siders distortions of the axisymmetric basis which is essen- 
tially our approach here, in a different context. A sample 
set of beam basis functions are shown in Fig. [2j 

We argue that this basis has several advantages for 
representing beams, and in particular we contrast it with 
Zernike polynomials, which are commonly used as a ba- 
sis set to represent distortions on the surface of an opti- 
cal element, and are sometimes used to represent beams. 
Elliptical Hermite-Gauss functions, by design, reproduce 
an elliptical Gaussian beam with the first basis coefficient, 
compared to Zernike polynomials, which take many more 
components to approximate it. Elliptical Gaussian beams 
are a frequently used and well-studied approximation for 
off-axis CMB instruments. Hermite-Gauss functions are or- 
thogonal over the whole plane while Zernike polynomials 
are limited to a disk, the reason they are convenient for 
distortions on circular apertures. This is not a problem for 
the Zernike basis, if we scale the disk larger than the area 
where we consider beam data, but the size of the disk we 
should use is ambiguous, particularly for beams with a wide 



variety of sizes. This same scale ambiguity is constrained in 
the our elliptic Hermite-Gauss basis by fitting the ground 
state to the beam. 

We checked that most of the simulated Planck beams, 
and the ones most important for CMB measurement, 
may be represented compactly in the elliptical Hermite- 
Gaussian basis (again, see Fig. [2j). However, the multi- 
moded horns (at 545 and especially 857 GHz) do not closely 
resemble elliptical Gaussians, and this basis is less effec- 
tive when a small number of basis functions are employed. 
Because the basis is complete, these beam can be repre- 
sented by pushing to higher order, but it may be inefficient 
and unappealing to do so. Possibly, a different basis may be 
more successfully applied to these channels, but we have not 
explored this. For the channels where this basis does work 
well, the beam data should always be tested with higher 
order modes than the optics simulations require, checking 
for evidence of unexpected features at very low signal. 

Elliptical Hermite-Gauss basis functions are written as 



$„,„, (X) = ^";lf;j;f"-/^^f exp (-X' . x72) 



(9) 



where H^ is the Her mite polynomial of order n. The pa- 
rameters of the best-fit elliptical Gaussian are encoded in 
the transformation 



X 



^. 
cr: 



1 R-i(V')(x-xo), 



(10) 



which offsets the beam position by xq, rotates the beam 
through tjj, and scales the ellipse axes by 



CTx = t ^/^ X 6>FWHM/\/81og2 
CTy = t^l^ X 6>FWHM/\/81og2, 



(11) 



where t is elliptical Gaussian's axis ratio and ^fwhm is the 
geometric mean full width at half maximum. Then we can 
expand the beam as 



5(x) =^5n<^n(x), 
n 

where we re-index the eigenmodes with 

(ni +n2)^ +ni + 3n2 



(12) 



(13) 



With the normalization used here, convenient for 
beams, the maximum of the ground state function is unity, 
and the orthogonality relation is 



/ d^X <l>^(x)<I>n(x) 



^^FWHM^^^ 



8 log 2 
So to recover the basis coefficients 5^, one integrates 



8 log 2 

^^FWHM 



/ 



£-x $„(x)B(x). 



(14) 



(15) 



Since Planck only samples the sky (Fig. [T]), these inte- 
grals must be approximated as sums over samples 



/d^.(...)«4E(-) 



(16) 



computed from the N detector samples from the area 
A surrounding the planet. In this sampling, the modes 
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Fig. 2. Basis suitable for decomposing elliptical beams. At left: the orthogonal and complete basis functions, where the 
ellipse marks the half-maximum curve Gaussian portion of the basis (bottom left). Every subpanel to the right or up 
increments the index on the Hermite polynomial perpendicular to or along the ellipse's major axis. At right: efficiently 
and inefficiently represented sample Planck beams, and their decomposition into the basis, limited to ni + n2 < 20. 



are only approximately orthogonal. The approximation be- 
comes poor if the typical size of the gap in the sampling is 
large compared to the FWHM of the Gaussian on which the 
functions are based, or compared to the scale of oscillations 
in the highest frequency modes considered. Uncorrected, 
this will lead to biased estimates for the beam coefficients. 
We address this by computing, to some maximum rele- 
vant mode, the symmetric overlap matrix of the basis func- 
tions on the sampled sky. 



h 



8 log 2 A 



^ ^m{^i)^n{'^i)' 



(17) 



For the plane subsampling, beams, and maximum mode 
for a typical Planck beam, this is a dense matrix, with 
diagonal entries with values near 1 and off-diagonal entries 
ranging from roughly —0.2 to +0.2. To get unbiased beam 
coefficients, we use the inverse of this overlap matrix to 
deconvolve, yielding an estimate of the basis coefficient: 



= E 



81og2 A 



y'$„(xi)di, 



(18) 



where the di are the time-ordered data, consisting of signal 
and noise, 



di 



z^n(Xi) 



(19) 



a 



O 



le+07 r 




15 20 
Beam FWHM 

Fig. 3. For varying re-pointing steps, the overlap matrix 
/ condition number for modes with ni + n2 < 10, for the 
Planck beams. 



noise in the beam coefficient reconstruction proportional to 
the condition number. Thus study of the overlap integral 
provides a means to evaluate quantitatively the effect of a 
given scanning strategy on the quality of the beam recon- 
struction. 



In the limit of a well sampled beam, the sums closely ap- 
proximate the continuous integrals, and Imn approaches the 
Kronecker Smn- Adding zero mean noise does not bias the 
estimated basis coefficients in the ensemble average. If the 
noise is white, characterized by a covariance Cov(ni,nj) = 
a'^Sij, then the covariance of the estimated basis coefficients 



Cov{Sm,Sn) 



81og2 A 2 1 

(J i^^^ 



'^^FWHM 



(20) 



In the linear operation which corrects for the sampling 
(represented by I~^ in Eq.[T8]), the condition number of the 
overlap matrix quantifies how well a particular sky sam- 
pling supports a set of basis functions. For a range of beam 
sizes and pointing realizations, we have plotted the over- 
lap matrix condition number in Fig. [3l Sparser sampling, 
for a fixed overall number of sampling points, boosts the 



2.8. Beam window functions 

For both our reconstruction techniques, our understand- 
ing of the beam impacts the cosmological interpretation of 
the CMB power spectrum: we deconvolve the beam to get 
an unbiased estimate. Define the beam window function as 
the function bf such that the ensemble average of the mea- 
sured power spectrum relates to the true power spectrum 
as (Ci) = bfCi. To make an unbiased estimate of the power 
spectrum, we divide the observed spectrum by bf at each 
multipole. 

For an asymmetric beam and an arbitrary scan strat- 
egy, computing the window functio n efficiently remains 
an open research q uestion (see e.g. iHinshaw et al.l 120071 : 
IShimon et al.ll2008[ ). However, we can make progress with 
suitable approximations. Over the course of a single survey, 
the Planck scanning strategy is approximately pole-to-pole. 
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so that over most of the sky (except near the poles) the 
beam strikes at a common orientation. The single orien- 
tation window function is readily computed in the flat sky 
approximation, where we can express the beam convolution 
as a multiplication in Fourier space. With 2-d wavenumber 
1, the measured signal is 5(1)5(1) and we write the mea- 
sured power spectrum as an average over azimuthal angle 
on the Fourier plane: 



^f (true) = 0.1 0.01 0.001 



Ci 



1 

2^ 



27r 



dl^B{\)s{\)B\\)s^{\) 



(21) 



Then the ensemble average yields the window function, 
{Ci) = ^fj dl^B{\)B*i\){si\)s*{\)) 



^jJdl^B{\)B*{\) 



Ci = bfCi, 



where the term in brackets gives window function. We en- 
force the normalization 6/ ^ 1 as / ^ 0. (In practice we 
normalize the lowest /-bin.) Compared to the true scanning 
strategy for Planck, this type of approximation imposes an 
error of fractio n of a percent on to the temperature window 
function (^Ashdown et al.l |2QQ9[ ) . If we include a detector 
transfer function in our fitting, we additionally convolve 
the beam map with a spatial filter representing the effect 
of the detector's response during a constant velocity scan, 
then deconvolve a filter based on the (slightly-different) fit- 
ted transfer function. 

The mistake we make when we deconvolve an approxi- 
mate beam is the ratio of the reconstructed window func- 
tion to the true window function: 



rf = bf/b. 



Ltrue* 



(22) 



The fractional error in the power spectrum is then 

ACi/Ci = rf-l. (23) 

Fig. m show this quantity for limited numbers of eigenfunc- 
tions of the non-parametric models in one of the 100 GHz 
channels (assuming no noise and dense sampling of the 
beam). Elliptical Gaussians fit to the simulated beam in 
real space make a poor approximation for computing the 
window function. 

For a set of Monte Carlo realizations of the fitted beam, 
we can compute the set of window function ratios, and char- 
acterize the covariance matrix due to beam reconstruction 
errorfi 



Cov(Q„QJ 






{rl,-l){rl,- 



(24) 



2 



!)• 



The covariance matrix is dense, strongest on the diago- 
nal and smoothly dropping away from it. Beam errors are 



^ One peculiarity of our approximation is that only three pa- 
rameters of the parametric model contribute to the window func- 
tion. It depends only on the scale (beam solid angle) and two 
components of shear. The other four parameters in the model 
do not contribute: the amplitude is an overall calibration; the 
position offsets translate to phases in Fourier space, and cancel 
in the multiplication of complex conjugates; and the rotation is 
integrated out. Indeed the beam covariance matrix in this case 
has only three substantial eigenvalues. 
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Fig. 4. Bias in the power spectrum caused by limiting the 
number of basis functions for beam reconstruction, showing 
the poor result from the elliptical Gaussian beam approxi- 
mation (m = n2 = 0). For HFI 100 GHz. 



highly correlated across multipoles. We find the correlation 
coefficients by normalizing the diagonal of the covariance 
matrix. For the parametric model at 100 GHz, the correla- 
tion coefficient between high (/ ~ 3000) and low (/ ~ 100) 
multipoles is very strong, at least 0.88. This is important 
for the measurement of cosmo logical parameters, such as 
Us or r, which couple across large multipole ranges. 

Some caution is appropriate regarding bias in the re- 
covery of the window function. The window function is 
a quadratic function of the beam. Therefore a procedure 
which reconstructs an unbiased estimate of the beam can 
produce a biased estimate of the window function. For our 
beam recovery techniques, we in practice detect little bias 
in the ensemble of window functions for channels at 100 
GHz and below. In higher bands, with smaller beams and 
higher signal-to-noise, the mean of the Monte Carlo ensem- 
ble (per /) slowly oscillates near, but slightly above or below 
the true window function. In principle a correction to this 
bias can be folded into the power spectrum analysis, but 
here we allow the bias to persist, to see if it has any effect 
in the cosmological analysis. 



3. Results 

Here we compile results for beam fitting and the subsequent 
errors imposed onto the CMB power spectrum. We take a 
single Jupiter crossing as our baseline case for beam fitting, 
and consider a case with destriped 1// noise, with no con- 
tribution from large-scale (/ < 250) CMB. (Confusion from 
signals on the sky can be removed because every region of 
the sky is re-observed at 7- month intervals). For each fre- 
quency band, we use one model LFI or HFI beam. Unless 
noted, we assume nonlinearities in the detector response 
have been corrected before processing. 

3.1. Parametric model 

The errors on the recovered model parameters are shown in 
Table[3]for 7- and 8-parameter beam models, based on 1280 
Monte Carlo steps with 2000 fitting iterations each (3000 in 
case of 353 GHz). The quality of the parameter recovery is 



Kevin M. Huffenberger et al.: Measuring Planck beams with planets 



Band 


<7cross 


<7co 


ctfwhm/FWHM cf^ 


^7+ 


^7x 


(J a/ A 


(Jt/t 


(GHz) 


(10-^ 


arcmin.) 


(10-^) (10-2 deg 


) (10-^) 


(10"') 


(10-^) 


(10-') 








Low-Frequency Instrument 








30 


218 


344 


937 267 


548 


1270 


164 




44 


125 


152 


572 53.1 


143 


118 


92.1 




70 


26.6 


35.6 


418 34.3 

High- Frequency Instrument 


38.9 
(known r) 


141 


51.6 




100 


7.24 


25.3 


157 13.0 


19.2 


36.1 


10.1 


0.00 


143 


3.74 


6.76 


69.7 4.30 


4.68 


6.17 


4.77 


0.00 


217 


0.473 


2.26 


23.4 3.86 


9.96 


3.86 


2.26 


0.00 


353 


8.96 


9.81 


13.3 4.72 


22.9 


2.86 


1.83 


0.00 


545 


0.295 


0.167 


5.15 0.167 


0.763 


0.521 


0.866 


0.00 


857 


0.351 


0.153 


1.87 0.287 
High- Frequency Instrument 


0.454 
(fitting r) 


0.239 


0.398 


0.00 


100 


7.13 


37.3 


144 13.8 


18.7 


37.7 


33.0 


96.3 


143 


4.41 


15.6 


80.8 5.36 


6.91 


7.10 


18.3 


88.8 


217 


0.494 


4.11 


30.0 5.03 


12.3 


4.75 


6.18 


26.7 


353 


49.3 


52.0 


72.2 25.8 


126 


20.7 


13.7 


34.1 


545 


0.300 


0.341 


5.49 0.172 


0.828 


0.521 


0.970 


6.25 


857 


0.351 


0.164 


1.87 0.287 


0.454 


0.238 


0.412 


1.38 



Table 3. Standard deviations on parameter estimates from a single Jupiter observation. The beam offset error in the 
cross-scan direction is across, and similarly for ctco- The beam angle ?/;, shears 7+ and 7x, and amplitude A are defined 
in Sec. 12.51 The fractional FWHM error is derived from the error on the scale parameter 5, also defined there. The time 
constant r is defined in Sec. 12.41 



exceptionally good, due to the very high signal-to-noise on 
Jupiter. Comparing channels, the relative quality of the fits 
is a complicated interaction between Jupiter's signal, the 
detector's noise, and the beam size (small beams concen- 
trate the signal but yield fewer useful data points). Except 
for 353 GHz, the quality of the fits tend to improve with 
increasing frequency band, largely due to the increase in 
Jupiter's signal at high frequencies. Repeated observations 
of the planets during subsequent surveys reduce the errors 
roughly as expected for independent observations. 

The errors at 353 GHz, although quite small, are puz- 
zlingly larger than the errors at 217 GHz and 545 GHz, es- 
pecially when fitting for a time constant. This seems to be 
due to the x^ minimization getting caught in local minima 
away from the global minimum, which creates a population 
of outliers in the Monte Carlo ensemble of fitted parame- 
ters, driving up the errors. These outliers represent ~ 5% of 
samples and are seen only in the simulated 353 GHz beams. 
Although the other beams span a large range of beam sizes, 
signal-to- noise, and time constant duration, none show any 
obvious population of outliers. As noted in Fig. O several 
353 GHz beams show this behavior. The outliers do not fall 
into any well-separated population which make them easy 
to cut, and we have not found a way to ehminate them 
robustly. 

The 353 GHz case (signal, noise, and time constant) run 
with the similarly-sized 217 GHz beam does not show a 
large population of outliers. Visual inspection yields noth- 
ing obviously wrong with the 353 GHz beams, which are 
formatted the same way as the other HFI beams, and the 
timelines the 353 GHz beams produce in our pipeline are 





4.22 




4.2 




4.18 




4.22 






a^ 


4.2 


1 


4.18 

4.22 

4.2 




4.18 




4.22 




4.2 




4.18 



-+- 



wu 



HIT 353-3a 



h- 



HFI 353-4a 
1 h- 



JHffWt>*YJL4^ 






HFI 353-5a 



-+- 



HFI 353 signal + noise, HFI 217-5a beam 



200 400 600 800 1000 
Independent MC realization 



1200 



Fig. 5. Parameter fits for the detector time constant for 
the HFI 353 GHz channel. Several 353 GHz horns show 
outliers (top three panels) not seen in the other channels. 
The outliers are not seen when a simulated beam from the 
217 GHz channel is substituted into the 353 GHz simulation 
(bottom panel). 



also unremarkable, so it remains unclear why these outliers 
occur. Even with the outliers, the beam fits for 353 are still 
quite good, only suffering by comparison to the spectacular 
results from the neighboring channels. 

For all channels, the corresponding errors on the power 
spectrum due to the window function uncertainty are de- 
picted in Fig.[6l Our ensemble of window functions presents 
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Fig. 6. Errors on the window function using the para- 
metric beam model. At each multipole, 68% of the fitted 
Monte Carlo window functions recover spectra closer to the 
true power spectrum than the indicated line. Lines are cut 
off where the window function falls to 1%. For the HFI 
(bolometer) channels thinner lines of the same color and 
type denote fixing the time constant before fitting, showing 
smaller errors. 
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Fig. 7. Errors on the window function, like Fig.[6l but using 
the non-parametric beam model limited to ni + n2 < 20. 



Because we have some knowledge of the beams from 
ground tests and numerical models, these parametric and 
non-parametric cases should bracket the range of reason- 
able possibilities. 



a slightly biased estimate of the true window function (sec- 
tion 12. 8p , so instead of a standard deviation, we plot a 
contour which bounds the error for 68% of the window 
functions in our ensemble, including both bias and dis- 
persion. Errors are strongly correlated between multipoles. 
Depending on multipole, prior knowledge of the detector 
time constant improves the errors on HFI channels 100- 
353 GHz by a factor up to 2-3. The errors on the higher 
frequency channels are less affected. In general, the recov- 
ery of the window function for the higher frequency bands 
is exquisite, but this is a result of the rigidity of the model, 
due to the small number of parameters. 



3.2. Non-parametric decomposition 

The non-parametric model (Fig. [7|) has notably larger er- 
rors. This is due to the flexibility of this model compared to 
the parametric one. For computational efficiency, we limit 
the basis coefficients to those with ni + n2 < 20. At that 
refinement, the 857 GHz channel in particular is poorly re- 
solved by this basis, leading to large errors in the window 
function. 

If the parametric model, where only a handful of num- 
bers can describe the beam, describes our best case realis- 
tic scenario, the non-parametric model, requiring no prior 
knowledge of the beam, represents our most conservative 
case. For the bands 100-217 GHz, which are the most im- 
portant in terms of raw sensitivity to the CMB, one of the 
twice-annual crossings of Jupiter should yield a measure- 
ment of the beam window function to within 0.3%. This 
compares to r oughly 0.5% from 5 years of Jupiter mapping 
from WMAP (|Hill et al.l [20091 ). Over the mission lifetime, 
Jupiter will be visible about four times. Errors will decrease 
somewhat faster than the square root of the number of ob- 
servations in the non-parametric model, because filling in 
the plane constrains which modes can contribute, improv- 
ing the overlap matrix. 



3.3. Detector non-linearity 

For the parametric model, we evaluate an extremely con- 
servative method for dealing with the nonlinear gain of the 
HFI, excluding data where the gain deviates from linear- 
ity by more than the rms noise per sample. In practice, the 
HFI analysis pipeline will correct the data for the nonlinear 
response. 

For most channels, this cut removes the peak region of 
Jupiter and Saturn, but keeps the central region of Mars, 
which provides information on the beam's peak, so we in- 
clude all three in the fit. Separate amplitudes are fit for 
each of the planets, and are effectively marginalized out 
in the computation of the window function. Only the 353 
GHz channel has significant nonlinear response at the peak 
of Mars; for this channel only we additionally include ob- 
servations of Uranus and Neptune to aid the fitting. The 
corresponding error in the window function is shown in fig- 
ure [H based on 384 Monte Carlo simulations per channel. 

Although we are including more planets, the exclusion 
of the peak of the beam on Jupiter, where signal-to-noise is 
highest, boosts the noise significantly. In this case the error 
on the window function on the small scale end of the beam 
is raised by a factor of 2-25, with 545 GHz the least and 
217 GHz the most affected. Correcting for the non-linearity 
will allow much better performance. 

3.4. Noise and destriping 

Destriping to reduce noise is an important step in beam re- 
construction. We ran cases with unfiltered low frequency 
noise, which manifests as stripes across the face of the 
beam, and this will increase the noise in the beam fit 
or decomposition. The non-parametric model is somewhat 
more sensitive to low-frequency noise than the paramet- 
ric model, because the correlated noise will be mistaken 
for the true structure of the beam. Larger beams are also 
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Fig. 8. Errors on the window function, like Fig.[6l using the 
parametric beam model, but including the nonlinear gain 
of HFI and multiple planets. 



more affected than small beams, because they take longer to 
transit the planet. For both beam reconstruction models, 
low frequency noise increases the window function errors 
at 30 GHz (FWHM 32') by nearly two orders of magni- 
tude and at 217 GHz (FWHM 6.5') by less than an or- 
der of magnitude. The parametric model errors increase at 
857 GHz (FWHM 4') by several tens of percent, but the 
non-parametric model's errors at 857 GHz are driven by 
the poor decomposition into basis coefficients, and are not 
much affected by the destriping. 



4. Implications for cosmology 

We explore the tilt of the scalar perturbation spectrum with 
the likelihood in a simplified case where the likelihood func- 
tion is analytic. We take slices through the likelihood, mod- 
ified by the beam errors, and leave a fuller Markov Chain 
Monte Carlo evalu ation of the likelihood to other work (see 
iRochaet al"] |2QlQ^ 

For a theoretical power spectrum, Ci = l{l -\- l)C^/27r, 
given a beam deconvolved data spectrum Vi , which includes 
isotropic noi se with power spe ctrum A/^, the full-sky likeli- 
hood is (e.g. iBond et al.ll2QQQ[ ) 



2\og C{Vi\Ci) = 

^(2/ + l) 



(25) 



\og{Ci^M)^Vi/{Ci^M) 



We assume a flat prior on Ci , so that the posterior probabil- 
ity is proportional to the likelihood: P{Ci\Vi) ex C{Vi\Ci). 
For the beam deconvolved noise spectrum, we take 



M = 



1(1 + 1) Airal 



2nbf 



tsurv 



(26) 



which is exact for uniform white noise, and depends only 
on the time domain noise variance (<j|, see Table [2]) and 
the duration of the survey (Tgurv) which we take as one 
year. To check the likeHhood for individual detectors, we 
construct noise power spectra for one detector at a time 
(Figure [9]). 
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Fig. 9. Noise power spectra compared to the CMB temper- 
ature power spectrum, for single horns in each of the Planck 
channels, assuming uniform white noise (from Table [2j) and 
1 year survey duration. (Noise figures for the full Planck 
focal plane, with 74 detectors, will be much lower.) 



4.1. Gaussian beam model 

To gain intuition on the magnitude of the beam impact 
on parameters, we use a 1-parameter symmetric Gaussian 
model for the beam window function, bf = exp(— constant x 
/^), where the constant describes the width of the beam. 
Deconvolving a mismatched beam yields a ratio of window 
functions which may be parameterized by the fractional 
error in the FWHM, which we denote A. For small errors, 
the window function ratio in this case is 



Ti =exp 



r 

2 A log a X — 



(27) 



which corresponds to rf = 1 + 2 A log a at the multipole 
defined hy bf = a, where the beam window function has 
fallen by a factor a. We include data only for / < /o.oi, that 
is, scales larger than where the beam window function has 
fallen to 1%. 

We fix all parameters except n^, which we expect to be 
the most sensitive to errors in the beam. For an incorrect 
window function, the data and the likelihood are distorted, 
as depicted in Figure [10] for the beam and noise of a 100 
GHz detector. We can quantify the bias in the likelihood 
by computing the distorted mean. 



,(A) = I dnsnsC{Vi{A)\Ci{ns)), 



(28) 



plotting the distance from the true value as a function of 
A, normalized by the error in n^ along that slice, given by 



/ 



dus {us 



ris^truef C(Vi\Ci{ns)). 



(29) 



Figure [TT] summarizes the likelihood bias for a detector in 
each of the channels, when the fitted beam is too small. 
The highest frequency channels (545 and 857 GHz) show 
no bias simply because the errors are so large, and are not 
displayed. The pivot point for the family of spectra in our 
slice is / ^ 570, so the channels (30 and 44 GHz) with larger 
beams, which weight lower / more strongly, show a negative 
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Fig. 10. Distortions to the likelihood slice for the scalar 
index rig when the fit for the FWHM of a Gaussian beam 
is too small by 0.1 and 1 percent, for HFI 100 GHz noise, 
with all other cosmological parameters fixed. 



using the window function ratio ri from Eq. [22j Assuming 
a flat prior on the window function, we may approximate 
the posterior probability distribution of window functions 
with our Monte Carlo ensemble of beam window functions. 

We find the fidelity of the reconstruction for the para- 
metric model is very impressive. In the bands most im- 
portant for the CMB, marginalizing over the ensemble in- 
creases the standard deviation of the likelihood slice for n^ 
by 6% at 100 GHz, 1% at 143 GHz, and leaves the errors at 
217 GHz essentially unchanged. We can also examine how 
the peak of the likelihood slice is shifted around for par- 
ticular realizations in the beam-fitting ensemble. For the 
parametric model, the rms peak shift is 0.6a at 100 GHz, 
0.3a at 143 GHz, and 0.1a at 217 GHz. 

The impact on the errors is more substantial in the non- 
parametric beam decomposition. The width of the likeli- 
hood slice for n^ is increased by 11% at 100 GHz; 9% at 
143 GHz; and 60% at 217 GHz. The ensemble rms bias in 
the peak of the likelihood is 1.2a at 100 GHz, 1.1a at 143 
GHz, and 1.7a at 217 GHz, indicating that the beam error 
is significant. 
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Fig. 11. Bias in the peak of the likelihood slice for n^ in 
units of the error, for single horns in a Gaussian beam 
model, when the fitted FWHM is too small. Horns at 545 
and 857 GHz are not shown because the noise is too large 
to effectively measure ng. 



bias in n^, while the others show a positive bias. For the 
three channels with the best noise, the beam requirement 
is strictest, and the fidelity in the beam required is striking. 
For example, at 100 GHz for a single detector, to limit the 
distortion in the likelihood slice to 0.1a, the FWHM must 
be known to 0.1% for a symmetric Gaussian. Note that 
our parametric fits (Table [3]) are achieving this precision 
in all nine channels. For a general beam, this means that 
the beam window function bf must be known to almost 
0.04% where it has fallen to 1 percent, and better at lower 
multipoles. 



4.2. Realistic beams 

We can marginalize the likelihood over the errors in the 
window function. 



/ 



d{ri} C{riVi\Ci{ns)) P{ri\ planet observation ), (30) 



5. Conclusions 

We have examined the problem of fitting beams for Planck 
to planet observations. Using a simple, but rigid, paramet- 
ric model, and a very flexible non-parametric model for the 
beam, we predict errors in the beam reconstruction from 
the focal-plane transit data with Monte Carlo simulations. 
As part of the development of the non-parametric beam 
decomposition, we showed how to evaluate the impact of 
a given scan strategy on the quality of beam reconstruc- 
tion. We note that elliptical Gaussian approximations fit 
to the simulated beams in real space produce substantial 
errors in the window functions, and should not be used for 
cosmological analysis. 

The errors are much smaller in the parametric model, 
but it is a toy model holding the place for a detailed opti- 
cal reconstruction of the telescope, which would probably 
require more parameters and provide less fidelity. The non- 
parametric model depends only on the data from planet 
scans, and requires no modeling of the telescope optics. 
Taking this as a pessimistic scenario for beam uncertainty, 
we project that a single transit of Jupiter should constrain 
the beam window function in the key 100-217 GHz CMB 
channels to 0.3%. This level of beam errors, however, will 
be a significant systematic for the measurement of the 
scalar spectral perturbation index Us as determined by a 
slice through the cosmological parameter likelihood. Other 
sources of uncertainty in the planet measurement, such as a 
calibration error or low frequency structure in the detector 
response, will raise the uncertainty in the final cosmological 
measurement. 

In this analysis, we have for simplicity excluded sev- 
eral effects which may prove important for the correct re- 
construction of the beam. These include uncertainties in 
the solution to the telescope's pointing, uncertainty in the 
planet's microwave frequency spectrum, relating to the dif- 
fering sizes of the beams for planets and for the CMB, time 
variability in the planet signal (including the Galilean satel- 
lites of Jupiter, which by themselves will be high signal-to- 
noise signals for Planck)^ and the finite size of the planet's 
disk. 
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